function [va vvr]=RAICAR1D_var(dire,n1,n2,n3)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%R.Magalhaes, Paris, July 2019
%Function to calculate the component order variance across the number of
%permutations (va) used as well as return a 4D matrix with the averaged matched
%components (vvr)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Must choose the range of values to build correlation matrix
%and step size (e.g between 5 and 400 permutations in steps of 5
%Usage: RAICARD1D_fin('directory',min number of permutations, size step,max number of permutations);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Loaded file can be created with matlab function RAICAR1D made available
%within the same set of functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This software is provided under a Free for Non-Commercial Use Only Agreement and a Non-Commercial use agreement
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%contact: ricardo.magalhaes@cea.fr

c=1;
load(strcat(dire,'/CRCM_mat1D400'),'set','sizes','index','crcm');
load(strcat(dire,'/CRCM_data400'),'data');
setall=set;sizesall=sizes;indexall=index;crcmall=crcm;
clear set, sizes, index, crcm;
for i=n1:n2:n3
i
set=setall(1:i,1:i);
sizes=sizesall(1:i);
index=indexall(:,1:i);
crcm=crcmall(1:sum(sizes),1:sum(sizes));
load(strcat(dire,'/RAICAR1Dresultsnew',int2str(i)),'comp');
vr=[];
d=[];
for ij=1:i
for j=1:size(set(1,ij).mat,1)
d=[d set(1,ij).mat(j,:)];
end
end
hd=hist(abs(d),200);
hdlog=log(hd);
hdsm=smooth(hdlog,125,'sgolay');
[m, in]=min(hdsm);

for j=1:length(comp)
j
ck=1;
if(comp(1,j).r(1,3)>in/200)
vv=data(1,1).file(:,:,:,comp(1,j).r(1,1));
vr(1,:,:,:)=vv;
v(1,:)=vv(:)';
ck=ck+1;
end
if(comp(1,j).r(2,3)>in/200)
vv=data(1,comp(1,j).r(1,4)).file(:,:,:,comp(1,j).r(1,2));
vr(2,:,:,:)=vv;
v(2,:)=vv(:)';
ck=ck+1;

end

for k=2:size(comp(1,j).r,1)
if(comp(1,j).r(k,3)>in/200)
vv=data(1,comp(1,j).r(k,4)).file(:,:,:,comp(1,j).r(k,2));
v(ck,:)=vv(:)';
vr(ck,:,:,:)=vv;
ck=ck+1;
end
end
if(isempty(vr));vr=zeros(1,48,48,39);end
vvr(j,:,:,:)=mean(vr,1);
vr=[];
va1=var(v,1);
va1(isnan(va1))=0;
v=[];
va2(j)=mean(va1);

end
vv2(c,:)=va2;
va(c)=mean(va2);
c=c+1;
end

end